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1 Introduction 



Hybrid Monte Carlo (HMC) |T| remains the most popular algorithm for simulation of Quan- 
tum Chromodynamics (QCD) including the dynamical effects of fermions. It is therefore 
imperative that we have a detailed understanding of how the computational cost of HMC 
varies as we change simulation parameters such as the lattice volume and the fermion mass. 
In this paper we present detailed analytical results for the computational cost of the gener- 
alised HMC algorithm for free field theory. We expect many of the results obtained to also be 
applicable to more general field theories where most of the "modes" are weakly interacting, 
for example the high momentum modes of asymptotically free theories such as QCD. 

In this paper we give detailed descriptions of the techniques we have developed to enable 
us to calculate fictitious-time evolution operators, Metropolis acceptance probabilities, and 
autocorrelation functions for arbitrary polynomial operators in GHMC simulations of free 
field theory. This allows us to minimise the computational cost of such simulations, and to 
compare the efficiency of various popular limits of GHMC. This paper brings together and 
extends many results that have been presented previously p|, |^, H, |^, H, |^, || . The techniques 



developed for this paper have also been used in several other papers |T0| |TT], |12| . 

The structure of this paper is as follows: in Section |^ we describe the Generalised Hybrid 
Monte Carlo (GHMC) Algorithm and discuss various limiting cases. Section ^ demonstrates 
that GHMC computations for free field theory are equivalent to GHMC computations for 
a set of uncoupled harmonic operators with judiciously chosen frequencies. In Section ^ we 
describe the generalised leapfrog discretisation schemes for the classical equations of motion 
in fictitious time, and introduce specific parameterisations of the evolution operators in or- 
der to facilitate the calculation of the Metropolis acceptance probability in Section ^. In 
Section ^ we introduce general techniques for calculating Laplace transforms of autocorre- 
lation functions of polynomial operators in free field theory, and present explicit results for 
arbitrary linear and quadratic operators for several of the limiting cases of GHMC, both 
for fixed- and exponentially-distributed GHMC trajectory lengths. Section |^ presents a de- 
tailed analysis of the cost of HMC simulations in the approximation that the Metropolis 
acceptance rate is unity, whilst Section ^ analyses the relative cost of GHMC, HMC, and the 
second order Langevin algorithm (L2MC) for non-trivial acceptance rates under some very 
weak assumptions. Finally Section |] discusses the somewhat related topic of how to opti- 
mise relative frequencies of updates and measurements in general Monte Carlo simulations. 
Some concluding remarks are contained in Section Several detailed technical results are 
relegated to Appendices. 



2 Generalised Hybrid Monte Carlo Algorithm 

For simplicity, we shall describe the Generalised Hybrid Monte Carlo (GHMC) algorithm for 
a theory of scalar fields, denoted generically by 0, with action S'(0). 
In order to generate the distribution 
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we introduce a set of fictitious momenta vr and generate the fields and vr according to the 
distribution 

le-^(^-)[#]H (1) 

where 

ff(0,vr) = l^vr,2 + 5(0), 

X 

and ignore vr. 

We begin by recalhng [|I3] that a Markov Process will converge to some distribution 
of configurations if it is constructed out of update steps each of which has the desired 
distribution as a fixed point, and which taken together are ergodic. The generalised HMC 
algorithm 0, is constructed out of two such steps. 

2.1 Molecular Dynamics Monte Carlo 

This in turn consists of three parts: 

i. Molecular Dynamics (MD): an approximate integration of Hamilton's equations on 
phase space which is exactly area-preserving and reversible; this is a map U{t) : 
(0, tt) ^ (</.', tt') with det = 1 and f/(r) = U-\-t). 

ii. A momentum flip F : tt — tt. 

iii. Monte Carlo (MC): a Metropolis accept/reject test 




F ■ U{(j), tt) with probability min 
{(f), tt) otherwise. 



This may be implemented by generating a uniformly distributed random number y G 
[0, 1], so that 

^ dye(e-^^ -y)= min(l, e""^^). 
The composition of these gives the MDMC update step 

[t') = [F- Uir) e{e-^^ -y) + ie{y- e~^% ( ^ ) . (2) 
This satisfies detailed balance because (F ■ f/)^ = I. 

2.2 Partial Momentum Refreshment 

This mixes the Gaussian-distributed momenta vr with Gaussian noise ^ 

/ tt' \ _ / cos d mid \ p / tt 
V ^ / " V -sin^9 cos^9 / \ ^ 
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If IT and ^ are Gaussian distributed, PG{TT)[diT] oc e ^^/^[rfyr], then so is Pg{^')'- 

poo roo 

Pc(7r') = / dU c/7rPG(O^G(7r)5(7r' + 7rcos^9-esm^?). 

The extra momentum flip F is included so that the trajectory is reversed upon an MC 
rejection instead of on an acceptance. 

We may combine the MDMC update of equation (|^) with the partial momentum refresh- 
ment of equation (j^) to obtain the following expression for their combined effect 



TT' 



/ 1 

COS sin -t} 
\ — sin -i? cos 



U{T)e{e 



-SH 



y)+F9{y-e 



-5H-S 



( <t> 
TT 



(4) 



2.3 Special Cases of GHMC 

Several well-known algorithms are special cases of GHMC: 

• The usual Hybrid Monte Carlo (HMC) algorithm is the special case where i} = 7c/2, i.e., 
the fictitious momenta p are replaced by the Gaussian distributed ^, the old momenta 
being discarded completely. The momentum flips may be ignored in this case as long 
as MCMD and momentum refreshment steps are strictly alternated. 

• 7? = corresponds to pure MDMC, which is an exact version of the MD or micro- 
canonical algorithm. It is in general non-ergodic. 

• The second order Langevin Monte Carlo (L2MC) algorithm of Horowitz |1l4| , |T5| cor- 



responds to choosing arbitrary but MDMC trajectories of a single leapfrog integra- 
tion step, r = 6t. This is an exact version of the second order Langevin algorithm 



17^, 18, 19 



The Langevin Monte Carlo (LMC) [|20[ algorithm has -d = tt/2 and r = 5t. This is an 
exact version of the Langevin algorithm. 



2.4 Tunable Parameters 

The GHMC algorithm has three free parameters, the trajectory length r the momentum 
mixing angle and the integration step size 5t. These may be chosen arbitrarily without 
affecting the validity of the method, except for some special values for which the algorithm 
ceases to be ergodic. We may adjust these parameters to minimise the cost of a Monte Carlo 
computation, and the main goal of this work is to carry out this optimisation procedure for 
free field theory. 

Horowitz pointed out that the L2MC algorithm has the advantage of having a higher 
acceptance rate than HMC for a given step size, but he did not take in to account that it also 
requires a higher acceptance rate to get the same autocorrelations because the trajectory is 
reversed at each MC rejection. It is not obvious a priori which of these effects dominates. 

The parameters r and d may be chosen independently from some distributions Pr{t) and 
Pm{^) for each trajectory (but of course they cannot be chosen in a way which is correlated 
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with the starting point in phase space). In the following we shall consider various choices for 
the momentum refreshment distribution Pm-, but we shall always take a fixed value for 
Phii'd) = SI"}) — i)o). The generalisation of our results to the case of other distributions of i) 
values is trivial, but it is not immediately obvious that it is useful. We choose each trajectory 
r length independently from some distribution Pr{t), as this avoids the lack of ergodicity 
caused by choosing a fixed trajectory length which is a rational multiple of the period of any 
mode of the system [^. This is a disease of free field theory which in interacting models is 
removed to some extent by mode coupling. 

3 Free Field Theory 

3.1 Complex Fields on a Finite Lattice 

Consider a d dimensional free scalar field theory on a \^ = L"^ lattice for which the expectation 
of an arbitrary operator Q{(f)) is defined by the functional integral 

(fi) ^ i I e-'m{^) 

where Z is chosen such that (1) = 1. For a complex field : — C the action is 

where d'^ is the lattice Laplacian 

d 

For the Generalised Hybrid Monte Carlo (GHMC) algorithm discussed in section ^ we in- 
troduce a set of fictitious momentum fields vr, and the corresponding Hamiltonian 

i7(0,7r)^| E 7rX + ^(0), 
x&Zl 

and we evaluate the functional integral 

(fi)^l|[rf0][rf7r]e-^(^'-)fi(0). 

For theoretical analysis it is useful to diagonalise the Hamiltonian by Fourier trans- 
formation to "real" (as opposed to "fictitious") momentum space. From the identity of 
example ( |B.1.1|) we have 

0p = (-^0)p = L 
0x = (^"V)x = L 



-d/2 



-j — 2iTip-x/ L 



-d/2 ^ ^2-Kip-x/L 



(5) 
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and similarly for the fictitious momenta. We obtain 



where the frequencies are 



p ' 



(6) 



(7) 



The phase space configuration (0, vr) is generated by GHMC with probability density pro- 
portional to 



[dcf>][dn] 



oc e 



e-^(^-) det (g) det (|) [4]^^] 
[d^][dn], 



where the last relation follows because the Jacobian is a (real) constant. We thus see that free 
field theory corresponds to a set of harmonic oscillators with a specific choice of frequency 
spectrum. 

3.2 Infinite Lattices 

For an infinite lattice Z'^ (instead of Z|^) the corresponding momentum space is the torus Sf 
(instead of ), and equations (H), (H), and (^ get replaced by 



bp = {J^(t))p 



(2 



-tp-x 



if (0, tt) = i7(0, tt) = i L (^P^p + ^p'/'p'^P y ' 



and 



a; 



(p)2 = m2 + 4^sin2^. 



where we made use of examples ( |B.1.2| ) and ( |B.1.3| ). 



3.3 Real Fields 

For real (as opposed to complex) fields : Z^ — R we define 

0l°) = Im0„ 
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Since = (pjf'^ and 0^.^ = — ^jj"-* these two fields are independent degrees of free- 
dom only on a subset of the momentum space p G Z^. We may choose to define the real 
momentum field (h : Z f M as 



V2 



if P > — p 
if p < —p 



where the ordering is arbitrary {e.g., lexicographic). The Hamiltonian is then 



-2 E {^1 + 4^1 

peZf 



because = u^. 



3.4 Spectral Averages 



Results obtained for a general set of harmonic oscillators will be expressed in terms of "spec- 
tral averages," which may be evaluated explicitly for free field theory. We use the notation 
cr(/) to denote the spectral average of some function / : ^ M, (t(/) = L~'^J2p£Z'l ^ (^)' 
The finite lattice results may further be expanded about the infinite lattice results to 
obtain a large volume expansion by means of the Poisson resummation formula (^) of 
section ( ^^21) , 

d / 1 X p27r 

dpf, cos{p^kL) f{p). 



E 

keZ 



Example 3.4.1 Consider the massless spectral average in one dimension (where the volume 
V = L'^ = L) 




where all the higher order terms vanish for large V because Jq'^ dp sm"^" ^ cos pkV = for 
kV > a upon integration by parts. This result is to be compared with the exact answer given 
in equation ( |3^ of Appendix 0. 
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4 Leapfrog Evolution 
4.1 Lowest Order Leapfrog 

We wish to consider a system of V uncoupled harmonic oscillators {4>p} for p e Zy. The 
Hybrid Monte Carlo algorithm requires us to introduce a corresponding set of "fictitious" 
momenta {vTp}, and the dynamics on this "fictitious" phase space is described by the Hamil- 
tonian 



The classical trajectory through phase space must obey Hamilton's equations (using the 
symplectic 2-form X^ili d-Kp A 



^Yp) 

dH 



dH 



dn. 



-^l<pp- 



p 



The leapfrog equations are the simplest discretisation of these which are exactly reversible 
and area-preserving, 



(j)p{6T) 



vrp(O) + np{0) \5t 
0p(O) + 0p(i5r) 5t 



0p(O) + vrp(|5r) 5t 



This leapfrog integration scheme is a linear mapping on phase spaceQ 



r. 



0(r + 5t) 
-n{r + 5t) 



Uo{5t) 



where the matrix 



1 - Ut' 



0(r) 
7r(r) 

5t 



-6t + i(5r'^ 1 - i(5r^ 



(9) 



satisfies det f/o = 1, as it must because the mapping is area-preserving. This lowest-order 
leapfrog integration agrees with the exact Hamiltonian evolution up to errors of 0(5r^), 



^HSt 



where 



H 



TT 



dH d dH d 



dn 



dcj) dn 



TT 



1 
-1 



TT 



(10) 



is the generator of a translation through fictitious time, and Uo^i is some operator on phase 
space. The error must be an odd function of 6t because leapfrog integration is reversible to 
all orders in 6t. 



^For the rest of this section we shall consider only a single oscillator with = 1, as everything will be 
diagonal in p and the frequency dependence can be recovered by dimensional analysis. 
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4.2 Higher Order Leapfrog 

Following Campostrini et al. |j2^, ^ we can easily construct a higher-order leapfrog inte- 
gration schemes with errors of 0{5r') by defining a "wiggle" 



Clearly this is area-preserving and reversible because Uq is. Using equation (0), we find 
that 



Ui{St) = exp ( - — —H6t 



I + f/o,i (2 - a?) +0{5t' 



if we choose (Xi = \/2 then the coefficient of 5t^ vanishes, and we can arrange the step size 
to equal that of the lowest-order leapfrog scheme by taking ai = 2 — ai. The explicit form 
for the first-order wiggle Ui{6t) is 



1 _ 4+4^+3^ ^^5 



144 



This construction can be iterated by defining 



Un-l 



'5t 



Oir. 



(12) 



which again clearly is area-preserving and reversible. This may be shown to have errors of 
the form 



1 + U„i5t 



2n+3 



0{5t 



2n+5^ 



(13) 



by induction on n: Assume equation (|T^) holds Vn < A^, then from equations ([l^ ) and ( |T3|) 
we find that 



f/jv(5r) = exp(^-^/J5r 



'N 



'5t\ 



2N+1 



0{6t 



2N+3\ 



which gives us equation (^) for = upon setting aN = ^^"v^ and = 2 — a^. 



4.3 Parameterisation of Leapfrog Evolution Operators 

In order to calculate the explicit form for a Molecular Dynamics trajectory consisting of 
t/6t leapfrog steps it is useful to parameterise the leapfrog matrices UniSr) in the following 
way. Reversibility requires that for any leapfrog matrix f/„(— 5r) = f/„((5r)~^ and area- 
preservation requires that det UniSr) = 1, hence 
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For the lowest-order leapfrog matrix of equation (^) we observe that the diagonal elements are 
even functions of 5t and the off-diagonal elements are odd functions of 5r; this property also 
holds for all the higher-order leapfrog matrices defined by equation ([T2|). We may therefore 
parameterise [/„ in terms of two even functions 



l + ft„i5r2"+2 + 0((5r2"+4) 



as 



T 



sin[K„((5r) 6t] 
—pniSr) sin[K„(5r) 6t] cos[/t„((5r) 6r] 



cos[K„(5r) St] 



(14) 



(15) 



We may easily compute the leading terms in the Taylor expansions of these functions. For 
the lowest-order leapfrog scheme we obtain 



1 



Po = 1 - § 
for the Campostrini "wiggle" 



IH H 

24 640 



5t^ + 0{5t^) = 1 + 0.0416 St^ + 0{6t^ 



128 



Sr^ + OiSr^) = 1-0. 125(5x^ + (5r^); 



Pi 



1 



32 + 25v^ + 20v^ 



1440 

f4 + 3v^ + 2^ 

1 + ^ 

288 



6t^ + 0{St^) ^ 1 - 0.06614 + 0{St^ 



St^ + 0{5t^) ^ 1 + 0.03804 (Jr"^ + 0(5/ 



and for the second-order "wiggle" 



«2 



1 + 



6t^ + 0{6t^ 



P2 



1 - 



/ 2;(l321294+1039220\/2+831376v^) \ 
+ (1146320+901600 v^+721280 v^) 
+ (958174+753620 v^+602896 v^) 
+ (811769+638470 \/2+510776 v^) 
\ +(726418+571340 v^+457072v^) ^ J 

881798400 

/ 2; (391582+312244 v^+247752v^) \ 
+ (348752+277664 v^+219456 v^) 
+ (293134+233308 ^+184248 f/i) ^ 
+ (243521+194042 -§^+153684 f/i) v^S 
\ +(211570+168880 v^+134352^) / « , a. 

^ 5t^ + 0(6t^) 

125971200 



0.02169 5 



r 



1 -0.04573 + 0(5r^ 



Approximate values for k„ i and pn,i are listed in Table |T[ We note in passing that the 
magnitude of these leading non-trivial coefficients do not grow rapidly with increasing n. 
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n 




Pn,l 







0.0417 


-0.1250 


-1.0280 


1 


-0.0661 


0.0380 


-0.9945 


2 


0.0217 


-0.0457 


0.2861 


3 


-0.0204 


0.0038 


-0.7422 


4 


-0.0258 


0.0118 


1.4128 


5 


0.0437 


-0.0483 




6 


-0.0137 


0.0371 




7 


-0.0528 


0.1178 




8 


0.1545 


-0.0813 





Table 1: The coeffcients and pn,i of equation (|l^ for n up to 8. The hmiting values of 
\ogiQ{5Hn) / X for r ^ oo with m = are also given for the values of n appearing in Figure |I[ 

4.4 Time Evolution Operators 

The parameterisation given in equation (|15D facilitates the computation of the time evolution 
operator ?7„(r, (5r) = f/„((5r)'^/'^'^ for trajectories of length r where, in general, r ^ 5t. We 
assume as an induction hypothesis that 



Un{k5T,6T) = Un{6T) 



r N ( r 1 sin[K„(5r) kSr] 
cosK(5r) Mr] ^ p^^^;^ ^ 

-pniSr) sin[K„(5r) kSr] cos[K„((5r) kdr] 



from this it immediately follows that UniSryUniSr)^ = Un{5Ty^^ using simple trigonometric 
identities. Expressing the result in terms of the trajectory length rather than the number of 
integration steps, r = kdr, we obtain 



r /r N T sin[/€„(5r)rl 
cosK(5r)r] 

-pniSr) sm[Kn{Sr) r] cos[K„(5r) r] 



(16) 



This time evolution matrix may be expanded about the exact Hamiltonian evolution as a 
Taylor series in St, 



(17) 



where from equation (|TT] 



COS r sm r 
— sin T cos r 



and 



/ . . / — sinr COST \ / 1 

t/n,i(r) = -p„,ismr ^^^^ sin r + '^"'^M -1 



It is instructive to compare equations (^) and (|T7p. 
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5 Acceptance Rates 

In this section we describe the calculation of the acceptance rate for MDMC (and thus for 
GHMC) for a system of V uncoupled harmonic oscillators ^ . The method of calculation 
is the same for both Langevin Monte Carlo and Hybrid Monte Carlo, and is independent of 
whether one uses lowest order leapfrog or a higher order discretisation scheme; the various 
algorithms differ only in the explicit form of the time evolution matrix [/„(r, (5r) given in 
equation ([T6|) . 

The Hamiltonian of equation (^) is a quadratic form 

^ _ 1 / UJp(j)p Y ( ^p4>p \ 

so the change in "fictitious energy" over a trajectory is 

1 ^ / UJp(l)p \^ / UJp(l)p \ 



6H^ ^ h{<P{t), 7r(r)) - h{<P, vr) = ^ ^ ( 1 



where 5M„ = f/Jf/„ — I, and we have abbreviated 0p(O) and 7rp(0) to 0p and Tip. Inserting the 
Taylor expansion of equation (|T7|) and noting that {e^'^)'^e^'^ = 1, corresponding to exact 
energy conservation for St = 0, we find that 



I + f/n,i5r'"+2 ^ 0((5r'"+^) ^ I + Un,i 5r'"+' + 0{5t 



T r 



_2n+2 I /n/x^2n+4-i 



2?i+4\ 



The probability of the change in energy 6Hn having the value ^ when averaged over the 
equilibrium distribution of starting points on phase space is 

where as usual the "partition function" Z is just the normalisation constant required to 
ensure that / dC,PsHn{0 — ^- is convenient to choose an integral representation for the 
5-function as a superposition of plane waves, as then 



'^'^ Jdet{I + i'nSMn) J 2vr 



It will prove useful to observe that the exact area-preservation property of the time 
evolution operator ensures that det f/„ = 1, and thus that 

exp tr ln(I + (5M„) = det U^Un = det det f/„ = 1, 
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and hence trln(I + (5M„) = 0. This imphes that the quantity trln(I + irjSJ^n) niust vanish 
not only for ir] = but also for it] = 1. Expanding the logarithm, and making use of this 
fact, we find that 



trln(I + iriSMn) = tvil ~ (dM^^) 5t^''+^ + 0((5r*"+^). 



(20) 



We may perform an asymptotic expansion of the integral over rj using Laplace's method. 



First we recast equation (^Of) into a form where the V dependence is explicit 
tr ln(I + ^r^5M„) = zr^{l - (5M^ ^) x + (i^) , 



(21) 



where we have introduced the variable x = l^^r^""*"^, and the spectral average cr(/) = ytrf, 
which has a finite limit as V ^ oo. In equation ( pT| ) the correction terms are negligible 
if the only important contributions come from the regions where rj^ -C ^"'^\/V. Laplace's 
method shows that such contributions are exponentially suppressed. Using equation (|2T| ) in 
equation ([T9|) we obtain 



Completing the square gives 

1 



27r 



exp 



X 



4 



and as (SHn) = J PsHniO^ ~ \ ^ (^-^n,!) ^ this may also be written as |2^ 



47t{6H) 



: exp 



A{6H) 



where we have supressed the index n. 

The average Metropolis acceptance rate is now easily found. 



P = 



mm ( i, e 





5H 



in (l 

d^PsHiO 



oo 

2 roo 



die 



rf^P5^,(0min(l,e-«) 

d^PsniOe-^ 
1 



erfc 



■{6H') 



The explicit form for 



{6Hn) = ^a(5My x = 2xp^_ia((sino^r)V"+^) 



(22) 



^ p=i 



2, ,4n+4 
P ' 



(23) 
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where 



1 

32 



2p\^ = ^ = 0.03125 
for the lowest-order leapfrog integration scheme, and 

(40 + 32^ + 25 v^) 



2p?,i 



2pL 



41472 



0.002894, 



/ 2;(221124374726+175505620552 v^+139299694184^) ^ 
+ (192732972517+152971330676 ■§/2+121414378108 v^) 
+ (167855113064+133226114824 ■?^+105742133216 v^) 
+ (145990198426+115872242816 v^+91967883784 v^) 
y +(126961318646+100768983952 v^+ 79980474344 v^) J 
793437161472000 



0.004183 



for higher-order "wiggles." 



5.1 One Dimensional Free Field Theory 



For the case of free field theory we can compute the spectral averages using the explicit form 
for the frequency spectrum. Using the methods of appendix we find that the spectral 
average 



an 



+ - 



[3-3Jo(4r)+4J2(4r)- J4(4r) 
2-2Jo(4r)+4tJi(4r) + J2(4r 

h(2t2-l)Jo(4r)+3tJi(4r 



m 



m 



(24) 



up to terms which vanish as ^ ^ oo. The finite volume corrections are given explicitly in 
appendix ^ For the higher order integration schemes the corresponding spectral averages 
are also given in Appendix 

The values for the logarithm of {SH)/x are shown in Figure |I|. The limiting values for 
r — i> oo for the massless case (m = 0) are given in Table ^ for m > the corresponding 
quantities diverge. 



6 Autocorrelation Functions 
6.1 Simple Markov Processes 

Let (00, 4>2, ■ ■ ■ , 4>N-i) be a sequence of field configurations generated by an equilibrated 
ergodic Markov process, and let (p,{(f))'j denote the expectation value of some operator f2 for 
(f) distributed according to the fixed point distribution of this Markov process. For simplicity 
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Figure 1: The logarithm of {5Hn)/x, where x = V6t'^^~^'^, for one dimensional free field 
theory in an infinite volume and with a mass m = 0.01 is shown for various orders of 
Campostrini wiggles {n = corresponds to leapfrog integration). See Table on page |Tl| for 
the limiting values of the curves as r ^ oo with m = 0. 
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we shall assume that (Q) = in subsections |0| and |6.2| . We may define an unbiased estimator 
Q over the finite sequence of configurations by 



N 



so (n) = ^Ef=i(w.; 



(Q) = 0. The variance of this estimator is 



1 



A^-l N-1 



j=0 j'=0 



1 



N-1 N-2 N-1 

{ E {mf) + 2 E E {m')m 

3=0 j=0 j'=j+l 

( o N-1 \ 



1 



N 



1=1 



where 



CM 



(25) 



is the autocorrelation function for Vl. If the Markov process is ergodic, then for large ^, 

|C^n(^)|<ALx = e-^/^-, (26) 



where Amax is the second-largest eigenvalue of the Markov matrix and N^^^ is the exponential 
autocorrelation time. If 3> N„^„ then 



= {l + 2v4n} 



N 



N 



1 + 



1 + 

/ ' cxp 



' CXf 



(27) 



where = I]fc=iOn(^) is the integrated autocorrelation function for the operator VL. 

This result tells us that on average 1 + 2^4^ correlated measurements are needed to reduce 
the variance by the same amount as a single truly independent measurement. 



6.2 Hybrid Stochastic Processes 

Suppose that a sequence of configurations (0(to), 0(^i), 0(^2), • • • , 0(^Ar)) is generated from 
0(to) as follows: the configuration 0(tj) is generated from 0(tj_i) by choosing momenta 
'7r(tj_i) as described in Section ( p.2|) , and integrating Hamilton's equations for a time interval 
Tj = tj—tj_i, where each trajectory length Tj is chosen randomly from the distribution Pfi^Tj). 



The autocorrelation function Cq defined by equation ( pSD may be expressed in terms of the 
autocorrelation function 



;f)(</.(t,))o(0(to)) 

Cq{Ti, . . . ,Ti) = 
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by averaging it over the refresh distribution 



poo 

Cnii)= dn...dT,PRin)...PR{n)Cnin,...,n). (28) 

Jo 

The integrated autocorrelation function then becomes 

°° POO 

An = y2 dn... dn Pr{ti) . . . Pr{ti) Cniji, . . . , r^). 

If we wish to determine autocorrelations in terms of the total elapsed fictitious time t = 
te — to = Y.k=i of the sequence of trajectories we may introduce yet another autocorrelation 
function by 

Cn{t) =J2 dn...dn Pnin) . . . PRin) 6{t - J2 Cnin,..., r,), (29) 
and in terms of this we find that 

/■oo 

An = / dtCn{t). (30) 
Jo 

The function Cn{t) has the advantage of giving the autocorrelations as a function of MD 
time which is approximately proportional to computer time. 

If we make the reasonable assumption that the cost of the computation is proportional 
to the total fictitious (MD) time for which we have to integrate Hamilton's equations, and to 
the volume V of the lattice^, then the cost €. per independent configuration is proportional 
to (1 + 2An)Vf/6T with f denoting the average length of a trajectory. The meaning of 
"independent configuration" was discussed in section |6TI| , and depends on the particular 
operator Q under consideration. The optimal trajectory length is obtained by minimising 
the cost, that is by choosing f so as to satisfy 

^ = ^ l + 2A„ + 2f^=0 

ar (IT 

1- - 



6.3 Autocorrelation Functions for Polynomial Operators 

In order to carry out these calculations we make a simplifying assumption: 

Assumption 6.3.1 The acceptance probability min(l, e"'''^) for each trajectory may be re- 
placed by its value averaged over phase space P^cc = ^min(l, e""^^)^; we neglect correlations 
between successive trajectories. Including such correlations leads to seemingly intractable 
complications. It is not obvious that our assumption corresponds to any systematic approxi- 
mation except, of course, that it is valid when 6H = 0. 

^For free field theory the volume V is just the number of uncoupled harmonic oscillators. 
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The action of a generalised HMC update on the fields 0, their conjugate fictitious mo- 
menta TT, and the Gaussian noise ^ for the trajectory under consideration from Equation (1) 
is 

UJ(j)' \ I UJ(f) 

tt' = 2t(r, y, 5H) tt 

where the matrix 21 depends on the trajectory length r, the noise rotation angle d (which can 
be chosen independently for each trajectory if we so wished), the uniform random number y 
used in the Metropolis accept/reject test, and the value of 5H. 

We may ignore the corrections of non-leading order in 5t to the MD evolution operator 
because for any given value of (Pace) of 0(1) there is a corresponding value of x which is 
also 0(1), and thus 5t is of order V^~i/(^"+*^). These corrections therefore only contribute to 
the autocorrelations through the acceptance rate itself at leading order in the large volume 
expansion. 



From the leading order contribution (|T^) we obtain 



2t 



with e+ = e 



I cos CUT sin CUT \ ( ^ ^ ^ 

— sin ujT cos ^ cos ujt cos ^ sin ^ \0j^+\ — cos ^ sin ^ | 6. 
\ sin UJT sin ^ — cos ujt sin § cos ^ j \ sin ^ cos ^ 



± ( e — y) , and uj being the appropriate frequency for each mode. 



6.3.1 Linear Operators 

We are interested in calculating the autocorrelation function for a general linear operator 

^] = 5]^]p0p (32) 



for a set of uncoupled harmonic oscillators {0p}. Such an operator is of course connected, 
meaning that = J2p^p{4>p) = 0- Its autocorrelation function may be expressed in terms 
of the autocorrelation functions for the individual harmonic oscillators themselves, for 

' - {n^) - ^^^^^^ 

and as the oscillators are uncoupled 

[(f)p{t)(f)p'{0)) = 6p^p>((j)p{t)(j)p{0)} + (1 - 5py)(0p)(0p') = 5py(0p(t)0p(O)), 

thus 
where 



{Mt)<ppm 
) 
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Let us proceed to calculate these single mode autocorrelation functions; while doing so we 
can drop the subscript p for notational simplicity. 

The average over the Gaussian distribution of gives (.^) = 0, so we can drop the last 
column of 2t, and since a new ^ is taken from a heatbath at the start of each trajectory we 
may drop the last row also. This leaves us with the basis = (a;0, vr), which is updated 
by = 2t*^^^$(^) where the matrix 



21(1) 



coscijr 



sm UJT 



sin UJT cos cos ur cos {} 



9, 



1 






cos'd 



6. 



(33) 



6.3.2 Quadratic Operators 

Let f2 be a generic quadratic operator for the set of uncoupled harmonic oscillators {0p}, 

= E ^pq(j)p(j)q (34) 

p>q 

whose connected part is l^c = ^ — (^) • The autocorrelation function for Qc is^ 



If we define the elementary autocorrelation functions for linear and quadratic modes to be 
Ctf)p{t) and C((^2) (t), then we may express C^X'^) in terms of them: 



p>q 



Higher degree polynomial operators may be treated in the same way. As before we calculate 
the single mode autocorrelation functions and again drop the subscript p while doing so. 

We cannot directly average the GHMC update matrix 21 over ^ as we did in the linear 
case, but we can linearise the problem by considering a homogeneous quadratic operator as a 
linear combination of the quadratic monomials ((a;0)^, Tiucj), vr^, ^uxp, ^vr, ^^). These quadratic 
monomials are updated by the symmetrised tensor product of the update matrix 21, 



/ MO' \ 



21®, 21 



7CUJ(f) 

V e 



^The integrated autocorrelation function for a disconnected operator diverges in general. 
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where the update matrix 21 ® ^ 21 


is exphcitly 








00 






221,^521^^ 


22l<^52l<^^ 










2l7r7r2t0^ 


2l.?2l^^+2t,^2l^^ 


2l7r?2l0^+2l7r7r2l^5 


or Of 








212 


22l.^2t^<A 


22l^52l7r7r 


212^ 










2t^^2t<^^+2l^^2t<^g 


21^521^^+215^21^5 










2l^^2l7r7r 


21^521^^+21^^21^^ 


215521.^+215^21^5 






I 


22l5.2l^<^ 


212 


2215^21^^ 


22I552I5, 


212^ 



/ 

As the update is now hnear we can average over the Gaussian distribution of ^ as before. 
The basis monomials become ((cj^)^, ttcjc/), vr^, 0, 0, {^'^) = 1), and this leads to three simpli- 
fications: 

• The fourth and fifth columns are multiplied by zero, and can thus be dropped. 

• The fourth and fifth rows are not of interest and may be dropped too, as ^ is chosen 
from a heathbath before the next trajectory, and we already know how the linear 
monomials (0, vr) are updated. 

• The last row may be replaced by (0, 0, 0, 0, 0, 1) as we know that C,' is Gaussian dis- 
tributed and thus (^'^) = 1. 



We are thus led to consider the update of the inhomogenous monomials of even degree in uj(f. 

\2 ,^ _2 1^ ^i^if^i^ jg updated by = 2l^2),|,(2) ^j^gj^g 



and IT alone, namely = ((tuf/))"^, ttcjc/), vr^, 1) 



21(2) 





91 2 




91 2 


9] 2 


\ 




2^7r</)2(-(^(^ 




2^7r7r2i07r 


21.521^5 






9l2 


221 2t ^ 


212 

TTW 


912 




V 











1 


/ 



COS CUT 

COS ur sin out cos -i? 
sin^ UT cos^ 




2 cos cur sin CUT 



sm UT 



(cos^ UJT — sin^ ujt) cos i} cos ujt sin ujt cos 



-2 cos cur sin cur cos^ ^ 




cos^ cur cos^ 







sin^ {} 
1 



/ 1 
















— cos?? 














cos^ 


sin^ 




V 








1 





(35) 



6.4 Laplace Transforms of Autocorrelation Functions 

We can extract more information about the autocorrelation function Co(t) than just the 
integrated autocorrelation function Aq. We shall discuss this further in Section ^. In order 
to do this it is convenient to compute the Laplace transform of the autocorrelation function 



FniP) = CCn{t)= / dtCn{t)e-^' 

Jo 

= J2 dn...dn PRin) . . . Pnin) Cn(ri, . . . 

Jo 



(36) 



r,)e-^"i . . . e-^"^ 
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We may generalise the results of sections |6.3.1j and |6.3.2| and observe that the update 
step for the vector of inhomogeneous monomials in ucp and vr, is of the form^ 

$W(r) =2l('^)(r)<l>('^)(0) 

The matrix 21*-'^-' is block upper triangular with with first block corresponding to the d + 1 
homogeneous monomials of degree d, the next block to the d — 1 homogeneous monomials 
of degree d — 2, and so forth. The dimension of the matrix, which is the number of even or 
odd inhomogeneous monomials of degree d or less, is i((i + 2)^ for d even and j{d + 3){d+ 1) 
for d odd. 

The connected operator {(f)'^)c = (p'^ — {(p'^) = ■ ^^'^^ where v = (l/uj^,0, ... ,0, {(p^)), 
and the ^-trajectory autocorrelation function is thus 

By virtue of assumption |6.3.1| we may replace each matrix 01^'^^ by its phase space average. 



.T 



C(0d)^(ri, . . .,Ti) 



■ V 



and the Laplace transform of the degree d single mode autocorrelation function is therefore 



V 



■ V 



■ V 



where 



2lW = dTPR{T)e-^^{Ql^''\T] 



(37) 



Summing the geometric series we obtain a simple expression for the Laplace transformed 
autocorrelation function. 



l + F( 



T 



■ V 



(38) 



In order to evaluate the expectation values, we need the Gaussian averages for the equi- 
librium probability distribution of equation (|l|) for the Hamiltonian given in equation (§): 



2k^ 



2k+l 







TC, 



2k 



and 

when p ^ q and r ^ s. 



7r2)(7rf 



After averaging over the distributions of and y which are chosen independently for each trajectory. 
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6.4.1 Linear Operators 

For linear operators the matrix of expectation values in equation ( ^8]) is 



$(1) $(1) 



1 
1 



The explicit form of the Laplace transformed update matrix is obtained from equa- 
tion (|33|) by first averaging over phase space, which replaces 6^ by P^^^iur) and 6_ by 
1 — P^^^(a;r), and then equation ( P7D gives 



^1,0 + ^0,0 ^0^1 

—Gq^i COS'S [Gifi — Gq q) cos-i? 



where 
and 



POO 

drPniT)e-^^il - P_(r))(cos^ry(sin^r)^ 



The Laplace transform of the linear single mode autocorrelation function is 



_ 2 



Gto + ^0,0 - [Gto - Go/ + Go^i cos ^ 



1 ~ ^1,0 ~ ^0,0 + ( ~'-^l,0 + '^0,0 + ^1,0 ^ C^Ofi + G'o,i ) cos^ 



which for the special case of = 7i/2 (HMC) simplifies to 

Gffi + ^0,0 



F^(/?;^ = 7r/2) 



1 ^1,0 ~ Gq q 



6.4.2 Quadratic Operators 

For quadratic operators the matrix of expectation values is 

/ 3 1 1 \ 
10 
10 3 1 

V 1 1 1 ; 



(J$(2) $(2)' 



and from equations (|35|) and (^7\ 



( ^tfl + Gqq 



where 



(39) 



(40) 





Gi^cos^ (C^o - ^0:2 - Co.o) cos^? Gticos^ 

2Gl^cos^^ (G^o + ^0,0) cos^ 7? Gcosin^T? 

^0,0 / 



cos2^9 















roo 

G,- = + GT^ = dr P^(r)e-^"(cos ut)' (sin cur)* 



(41) 
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The Laplace transform of the quadratic single mode autocorrelation function F(02)^(/3) is 



G 



2,0 



0,0 



'•-^2,0 ~ ^'-^1,1 T '-J0,2'-^0,0 T '-^0,2'-^2,0 "r '-^0,0 



cost? 



Go^2 + Go,o) (Go,o + + (cos^)2 
+ (^0,0 + G^o + ^0^2) (G^o' - + + 4G+i' - Go,o') (cos^)= 



1 ~ Gq — G 



2,0 



2,0 

Go,0 
2 



v_ 2 



0,2; 



"^2^0 + G'0,0 ~ 26*2^0^0,0 + ^0^2 "^0,0 + G'2^o)(cost9) 



/ ~ '-^0,2 ~ '-'O.O + <-'0,2*-^2,0 ^ '-^0,2'-^2,p \ 

+G^a2 G 



2,0 



'G'o,( 



V + 



'-^0,2'-J0,0 



'-^0,2"-^0,0 "T '-^2,0'-^0,0 '-^2,0 "T '-^2,0 



^0^2 ^0,0 



^0,0 + 2(70^2^2^0^0,0 / 



which for the special case of ^9 = 7r/2 (HMC) simplifies to 



(cos-i?)' 



(42) 



F(<^2)^(/5;^ = 7r/2) 



G 



2,0 



Go,( 



2,0 



0,0 



6.5 Exponentially Distributed Trajectory Lengths 

To proceed further we need to choose a specific form for the momentum refresh distribu- 
tion. In this section we will analyse the case of exponentially distributed trajectory lengths, 
Pr{j) = re~'^'^ where the parameter is just the inverse average trajectory length r = 1/r. In 
section |6.6| we will consider fixed length trajectories. 

We make the approximation that 6H and thus P^^ are independent of r (c.f., Figure |l|), 



so 



/^+ p n r^- _ p ■ 

^j,k acc^ j,kl j,k ~ V ^ acc)^ j,kl 



where Pace = -Pacc(^)- This approximation is only made in order to avoid unpleasant integrals 
which cannot be evaluated in closed form. The integral in equation (^) may be evaluated, 
and we find 



/^cxp 



EE 



J 



(i)i+'=(-l)-+L^^'JB 



^=-^0 .=0 VJ V) + 1)' + ^\3 + k-2^i- 2uf ■ 
where we have introduced the dimensionless quantity = uf and 



(43) 



B 



/?f + 1 



if iA; G N, 



6.5.1 Linear Operators 

The Laplace transform of the linear single mode autocorrelation function for exponentially 
distributed trajectories is obtained by substituting the explicit form for G^ into equa- 
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tion (0), and we obtain 



r(^/53+((l-2P_)cost9 + 3) r(3^ 
+ (2(1 - 2P_) co^d + (1 - P_)^2 ^ 3^ 
((1 - P_)V' + 1 - 2Pacc) cos?9 + (1 - P_)<^2 ^ ^3 



13^ + ((1 - 2P_) cost? + 3) r/?3 + (2(1 - 2P_) cos?9 + 3 + (^2) r^/^^ 

■ :i - P_)^2 ^ 1 _ 2P_) cost? + (1 + P_)</.2 ^ i\ ^3^ 



Pr(/3) 



+ (P_(l - P_)<^2 pQg^ ^ p^^^^2^ ^4 

Evaluating this at /? = gives the integrated autocorrelation function 

^..p _ ((1 - P.cc)V' + 1 - 2P..J cost? + (1 - P..J(^^ + 1 

^ P_(^2((l_p_)cOSt?+l) 

For HMC where t9 = 7r/2 we have 

''2 + 2r/3+((l-P_)<^2_^l)r2)r 



with the corresponding integrated autocorrelation function being 



P P 092' 

In the limit of unit acceptance rate, P^^^ = 1, we obtain 

(/3 + (1 — cos t9)r)r 



Fr(/3;P_ = 1) 



/32 + (l-cost9)r/5 + rV 
1 t9 

^r''(-Pacc = 1) = 5 • 

Finally, for HMC in the limit of unit acceptance rate 

(r + P)r 



P!-(/3;t9 = 7r/2,P_ = l) 



^2 + r/? + rV' 

and 

^7(t9 = 7r/2,P_ = l) = ^. 
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6.5.2 Quadratic Operators 



The Laplace transform of the quadratic single mode autocorrelation function for exponen- 
tially distributed trajectories is obtained by substituting the explicit form for G"^ into equa- 
tion and we obtain 



T-lCXp 



+ 



(3^ + [-{cos^f + (1 - 2P_) COST? + 4)r/5^ 

^ (-l + 2P_)(cos?9)3-3(_cost9)2 ^ 

^ +(3 - 6P_) COST? + (4 - 2P_)(^2 + Q 

(4P_ - 2)(cost9)3 + ((4P_ - 4)(^2 _ 3)(cost9)2 

+ (2(-l + P_)(P_ - 2)y.2 + 3 _ gp^^j ^ost9 I r'P 

+ (8-4P_)¥;2 + 4 

/ (_4(-l + P_)V-l + 2P_)(cosT?)3 ^^ 
+ ((4P_-4)(^2_^)(^Qg^)2 

+ (2(-l + P_)(P_ - 2)y;2 + 1 _ 2P_) COST? 
+ (4-2P_)y;2 + l 



{cos^f + (1 - 2P,,J COST? + A)r(3 



( (-1 + 2P_)(cost9)3 - 3(cosT?)2 
\ +(3-6P,,JcosT9 + 6 + 4¥p2 ^ 

/ (4P_ - 2)(cosT?)3 + ((2P_ - A)ip^ - 3)(cosT?)2 
+ ((4 - 4P_)^2 ^ 3 _ gp^^j pQg^ 

+ (8 + 2P_)y;2^4 
/ (-2(-l + P_)(P_ - 2)y.2 _ 1 + 2P_)(cost9)3 \ 
+ {-Aip^ -l){cos^f 
-2(P_ + 2)(-l + P_)^2 ^ 1 _ 2P_) COST? 
+ (4P_ + 4)v92 + 1 



V 

/ 2P_(-1 + P_)(^2(^Qg^)3 _ 2(C0ST9)2P_^2 

I -2P_(-1 + P_)(^2^ost9 + 2P_(^2 



Evaluating this at /5 = gives the integrated autocorrelation function 



/loxp 



4(1 - P_)V' + 1 - 2P_)(cost9)3 + (4(1 - P_)v^2 + iVposT?) 



+ 2(1 - P_)(P_ - 2)^^ - 1 + 2P\ COST? + 2(P_ - 2)v9^ - 1 



2P_(/^2(cost9 - 1)(cost9 + 1)((1 - P_) cost? + 1 
For HMC where t9 = 7r/2 these equations simplify to 

^2 + 2r/? + (2(2-P_)¥;2 + l)r2)r 



and 



A-.)Jt9 = 7r/2) 



2 P^rn 1 

+ 



P, 



2P.cc<^2 
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In the limit of unit acceptance rate, P^cc = 1, we obtain 



-^("^2) -Pace - 1) 



+ (- cos^9 + 2 - {cos^f)r(3 
+ ((cosi9)^ - {cosily - cos-d + 2(^2 _^ iy2 



and 



"^f^2> (-Pace — 1) 



/33 + (- COS ^9 + 2 - (cos 79)2)r/52 
+ ((cos^9)^ - cos^9 + 1 + 4v92 - (cos^9)2)r2/? 
+ (-2(cosi9)V^ + V)^^ 

(cos ■(9)3 - (cos'i?)^ - cos?? + 2v92 + 1 



2^2^cos?9)2 - 1) 
Finally, for HMC in the limit of unit acceptance rate 



F(7)^(/?;^ = 7r/2,P... = l) 



and 



A-/., (^9 = 7r/2,P_ = l) = l + 



1 



(45) 



2^' 



2 ■ 



6.6 Fixed Length Trajectories 



In this section we consider the case of fixed length trajectories, Pr(t) = 6{t — t). In this 
case we find that 

Cj^fe = -PaccCjr^fc and Gji^ = (1 — -Pacc)G'i,fc; 

without making any approximations, and from equation (^) we obtain 

G«:, = e-^^(cos¥.)^(sin¥,)'=. 

6.6.1 Linear Operators 

The Laplace transform of the linear single mode autocorrelation function for fixed length 
trajectories P^''(/5) is obtained by substituting the explicit form for into equation (^0]) , 

(-1 + 2P,,J cos^e-2^^ - (P,,, cosy^ + 1 - Pacc)e-^^ 

(1 - 2P_) cos^?e-2/3- + ((P_ + P_ cos if - 1) COST? + P_ cos(^ + 1 - P^^^)e-P- - 1 

Evaluating this at = gives the integrated autocorrelation function 

(1 - 2P,,J cos ^ + Pace cos V2 + 1 - -Pace 



For HMC where ■(9 = 7r/2 we have 



Pace(cos?9 + l)(cosv2 - 1) 



Pflx.^. n _ „/2^ _ (^aecCosy; + l-P.Je 

^ " / ^ " l-(PaceCOS^+l-P.Je-/3-' 
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with the corresponding integrated autocorrelation function being 
In the limit of unit acceptance rate, P^cc = 1, we obtain 



COS(y9e' 



COS 



and 



Afk (Pace — 1) 



^cos-i^cosf/) + cosv5)e — cost^e 2/3t _ ]^ • 
COS 'd — cos if 



cos ■(} COS V9 + COS if — COS ■{} — 1 

Finally, for HMC in the limit of unit acceptance rate 



and 



P«'=(/?;^ = 7r/2,P_ = l) = - 
A«^(^ = 7r/2,P_=1) = . 



cosy^e 



— cos 996 

cos ip 

1 — cos if 



-/3f ' 



(46) 



6.6.2 Quadratic Operators 

The Laplace transform of the quadratic single mode autocorrelation function for fixed length 
trajectories is obtained by substituting the explicit form for into equation (^2]), and we 
obtain 

(-P_(cos^)2-l + P_)e-^^ 
_e-3/5^(_i + 2P_) (cos 79)3 
+e-2/3^(-2P_ + 1 + 2P_(cos<^)2)(cos^)2 
+e-2/5^(P_ + P.cc(cos¥^)2 - 1) COST? 



TTlflx 



((-P...(C0S<^)2-1 + P_)(C0S^)2 

+ (-2P_(cos(/?)2 + 1) COST? - P_(cos(/?)2 - 1 + P_)e-^^ 
+ (^_e-2/3. + e-2/^^P_ + e-2/3^P...(cosy.)2 + e"3/5^ - 2e-3/^^P_) (cost?)= 
+e-2^"(-2P_ + 1 + 2P_ (cos 9^)2) (cos t9)2 
+e-2^^(P_ + P.ec(cosy^)2 - 1) COST? + 1 

Evaluating this at /? = gives the integrated autocorrelation function 

(1 - 2P_)(cost9)2 + 2cosT9Pa,,(cos<^)2 - P_(cos(^)2 - 1 + P_ 



(cost? — 1)(cost9 + l)(cosv9 — l)(cosv9 + l)Pa< 
For HMC where d = tt/2 these equations simplify to 

(P_(cosy.)2 + l-P_)e-^^ 



Pf;.) J/5; t9 = 7r/2) = - 



(P,ec(cOSV?)^ + 1 - Pacc)e" 



■I3f ■ 
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and 



P_(C0S(^)2 + 1-P,, 



P,,,(l - COS + cos(p) 
In the limit of unit acceptance rate, P^cc = 1, we obtain 



-Pace - 1) 



{cosifYe-'^^ + (cos^9)3e-3^^ 
-e~2^^(-l + 2(cos(^)2)(cos?9)2 



(cos'i9)^(cos(/3)^ + (—1 + 2(cosv9)^) cos-i? + (cosv9)^)e 



(cos(^)^ + e ^'^^)(cos'/?)^ 



and 



e-2/3r(_l + 2(C0S¥?)2)(C0S?9)2 - e^2/3r(pQgy,)2pQg^ _ ^ 

(cos-i?)^ — 2 cos -(9 (cos (y?)^ + (cosy^)^ 



'^""^ ' (cOS'i9)^(cOSV5)^ — (cOS^9)^ — (cOS(y9)^ + 1 

Finally, for HMC in the limit of unit acceptance rate 

(cOS(y9)^e~^^ 



and 



Pf;.)J/3;^9 = 7r/2,P_ = l) 



A^;2)J^ = 7r/2,P_ = l) 



1 — (cOS(y9)^e ^"^ ' 

(cos(^)^ 
1 — (cosy^)^ 



7 Autocorrelations and Costs for HMC with P^^^ ~ 1 

In this section we calculate autocorrelations for the special case ^9 = (HMC) and in the 
approximation where the acceptance rate is unity, P^^^ = 1. We shall consider more general 
cases in the following section. 

We begin with some general observations: 

• Integrated autocorrelation time. It is immediately obvious that the integrated auto- 
correlation time may be obtained from the Laplace transform (|36D by evaluating it at 

• Exponential autocorrelation time. For an ergodic Markov process we can write the 
typical asymptotic form of the autocorrelation function as 

Cn(t) ~ constant x e"*/*""^ for t ^ oo. 

This exponential autocorrelation time is a different quantity from A^^^p of Section |6.1| , 
which was defined in terms of Markov steps, t^xp can be extracted from Pn(/?) by 
considering its analytic structure in the complex (3 plane. Since t^xp governs the most 
slowly decaying exponential, it follows from the definition of the Laplace transform 
( p6D that Fq{P) will have its rightmost pole at Re/5 = — l/toxp- 
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• Dynamical critical exponent. One of the most relevant measures of the effectiveness of 
an algorithm for studying continuum physics is the exponent z relating the cost € to 
the correlation length ^ of the system, C oc 

• Optimal choice for d. For the GHMC algorithm we should minimise the cost by varying 
both r and d. For the case of quadratic operators with exponentially distributed 
trajectory lengths, for instance, the optimal choice of parameters when P^^ = 1 is to 
take — and r = ^d"^ — >■ 0. However, this ignores the fact that the cost does not 
decrease when we take r smaller than the 5t required to obtain a reasonable Metropolis 
acceptance rate. If we choose f^pt = 5t (L2MC) and the corresponding value for 'd^^^ 
we find that the cost is less than for the HMC case, but only by a constant factor. As 
the cost is only defined up to a implementation dependent constant factor anyhow we 
may conclude that generalised HMC does not appear to promise great improvements 
over HMC. The situation is more complex in the real world where P^^ ^ 1, which is 
explored in Section ||. 

7.1 Exponentially Distributed Trajectory Lengths 
7.1.1 Linear Operators 

The least uninteresting linear operator in free scalar field theory is the magnetisation M = 
X]x0(^)- From equations and (|^) this is simply the zero-momentum mode in Fourier 
space, uj = m, and is thus expected to evolve most slowly in fictitious time. 



Example 7.1.1 From equation ^3^), with Qpq = 6pQ, and equation ^4^), the Laplace trans- 
form of the autocorrelation function for the magnetisation is 



/3(r + /3) + m2 



As explained previously, the exponential autocorrelation time texp = — l/Re/3cxp corre- 
sponds to the rightmost pole in Fm, and this occurs when 



and therefore 



± 

exp 




2f 



1 / = if r < 

2^ if^>27^' 



^GXp 



where we have used the fact that r = 1/f. We observe that the exponential autocorrelation 
time is minimised when the average trajectory length is chosen to be f = l/2m. Note that 
M only couples to the zero momentum mode, and its relaxation rate determines t^^p in this 
case. 

The integrated autocorrelation function Am is given by 



Am = Fm(<)) 
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and we can minimise the cost of computing the magnetisation by making use of equation (j57[j. 
The optimal trajectory length is r^pt = \^/m, which corresponds to the minimum integrated 
autocorrelation function value Am {f opt) = 1/2. Note that the optimal trajectory length does 
not minimise the exponential autocorrelation time — they differ by a factor of 2\/2. 

The correlation length for this system is ^ oc 1/m, and our result indicates that the 
optimal trajectory length is f^p^ oc with a dynamical critical exponent z = 1. Indeed, if 
we were to choose an average trajectory length f oc then we would find that the cost per 
effectively independent configuration would grow as 



€ oc 



C," for a > 1, ,^ ^ cxD 
^2-° fora<l,^^oo. 



Keeping the average trajectory length fixed as the correlation length ^ increases, i.e., choosing 
a = 0, thus leads to a dynamical critical exponent z = 2. 

We can calculate the autocorrelation function CM{t) in closed form by inverting the 
Laplace transform. If we expand Fm in partial fractions 



A(l-A + 2r/?) A(l + A + 2r/?)' 
where A = yT^^"(2mfp, then it is easy to verify that the autocorrelation function is 
C,,{t) = J^{(l + A)e-(i-^)*/2.-(l-A)e-(i+^)*/2.| 

= ie-*/2^{cosh(^)+^sinh 



2f 



when < A < 1. When A'^ < we have 



C.,„^ie-fo.(f).ls.(f)), 

where A' = iA = ^ {2mfY — 1- We observe that there are oscillations in the autocorrelation 
function for long trajectories where f > l/2m. For the critical case A = we have 

Finally, the autocorrelation function expressed in terms of Markov steps is Cuin) = 
g-n/To^ yjliQj.Q iiiQ exponential autocorrelation time A^^xp = = l/ln[l + (mf)^]. 

7.1.2 Quadratic Operators 

Example 7.1.2 We obtain the Laplace transform of the connected autocorrelation functior^ 
for by setting Qpg = SpoSgo in equation ^^). From equation 

_ + 2r^(3 + r(3^ + 2m\ 

MllPJ ^3 + 2r/52 + (r2 + 4m2)/? + 2m2r' ^ ' 



^This is a synonym for the autocorrelation function for M^, the connected part of A'P 
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The integrated autocorrelation function is thus 



Am^=Fm2{0) = 1 + U— 

2 xniT 

Minimising the cost by means of equation ^3lD, we obtain f^pt = 1/ (^y/Srnj and thence 
^M2(^opt) =5/2. Again, the dynamical critical exponent is z = 1. The optimum trajectory 
length is different from that for M , which is to be expected. 

For exponentially distributed trajectory lengths the Laplace transform of the autocorrela- 
tion function is a rational function in (3 with the numerator of lower degree than the denomi- 
nator (see, eg, equations and ^jB^)), which implies that the autocorrelation function is a 
sum of exponentials. In this case the exponents are the roots of the cubic denominator, and 
they are either all real, or one is real and the other two are complex conjugates depending on 
the value of the mean trajectory length 1/r. This is shown explicitly in Appendix^. 



Example 7.1.3 Following equation ^B^), it is easy to write down the the Laplace transform 



of the connected autocorrelation function for the energy E = J2p h^f^"^ 



The integrated autocorrelation function for the energy is 

1^1 11 

^, = F,(o) = i + _i;-.i + _-.y) 



and the optimal trajectory length is fopt = \l o'l^'}^ /?>, leading to an integrated autocorrelation 
function value of AE^fopt) = 5/2. 

For one dimensional free field theory it remains to evaluate the spectral sum o'l^2 \ details 
of which are discussed in Appendix 10. In the physically interesting limit of small m and 
large V, we find cr^2 ~ Irn, and thus f^pt ~ l/yGm. Hence the dynamical critical exponent 
for the energy is i. 

For two dimensional free field theory we get z=0 up to logarithmic corrections (see Ap- 
pendix^. 



7.2 Fixed Length Trajectories 
7.2.1 Linear Operators 

From Section tlie Laplace transform of the autocorrelation function for the general linear 
operator of equation (^) is easily expressed in terms of equation p6|). The exponential 
autocorrelation time t^xp is related to the rightmost pole of the Laplace transform of the 
autocorrelation function, 

^cxp 1 / /3exp • 
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Equation ( ^B]) has poles in the complex (3 plane for 

/3 = — In coscUpf, 
r 

and hence ^ 

Re /3 = — In I cos cUpf I , 
r 

which leads to 

r 

~ max ■ 



exp 



In I cos uOpT I 



This diverges when Upf/ix G Z for any p, which just reflects the fact that the Hybrid Monte 
Carlo algorithm is not ergodic for these cases, as was first pointed out by Mackenzie [pT| . 
Perhaps this is most simply understood by considering the trajectory of the harmonic oscil- 
lator with frequency uOp in the {ujp(j)p, vr^) phase space. The Molecular Dynamics trajectories 
are circular arcs subtending an angle of Upf about the origin, and the momentum refresh- 
ment corresponds to a change of the Hp coordinate leaving the UpCpp coordinate unchanged. 
If ujpf is an integer multiple of vr the value of (pp can at most change sign, and thus this mode 
certainly cannot explore the whole of its phase space. 

Example 7.2.1 For the magnetisation the integrated autocorrelation function is 

cos mf 



Am = FuiO) 



1 — cos rriT 



If we minimise the cost we find that the optimal trajectory length corresponds to mf being an 
odd multiple of tc, and the worst case occurs when mf is an even multiple of n. Both cases 
just reflect the non- ergodic nature of the updating scheme discussed above: the fact that the 
optimal "ergodic" update occurs when mf is taken very close to an odd multiple of n is just 
an "accidental" consequence of the fact that (M) = 0. 

7.2.2 Quadratic Operators 

Example 7.2.2 In the case of the quadratic operator we find from equation (|6'. 6'. ^) that 



cos^ mr 



e'^^ — cos^ mr ' 



and thus Am^ = Fai^{0) = cot^ mf for fixed length trajectories. The optimal trajectory length 
TTtTa-pt ~ 1.3866. As is to be expected, the non-ergodicity at mfj-n G Z manifests itself as a 
divergence in Ami- 

Example 7.2.3 For the energy of one dimensional free field theory we obtain 



N p e^"^ — cos^ ujpT ' 



and therefore Ae = Fe{0) = jfJ2pCot'^ujpT which diverges whenever UpT / n G Z for any p. 
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8 Comparison of Costs for P^^^ ^ 1 



We wish to compare the performance of the HMC, L2MC and GHMC algorithms for one 
dimensional free field theory. To do this we shall compare the cost of generating a statisti- 
cally independent measurement of the magnetisation M and the magnetic susceptibility M^, 
choosing the optimal values for the angle d and the average trajectory length f. 

Following equation (^) we can minimise the cost with respect to d without having to 
specify the form of the refresh distribution. 

The next step is to minimise the cost with respect to the average trajectory length 
if = mf. Strictly speaking we should note that the acceptance probability P^cc is a function 
of f, but to a good approximation we may assume that P^cc depends only upon the integration 
step size 6t except in the case of very short trajectories, such as Langevin-type algorithms 
(see Figure |1|). The exact solution of the problem would clearly be highly transcendental. For 
this reason we shall find that although L2MC is a special case of GHMC and thus can never 
be cheaper, our "optimal" solution^ for very short trajectories (acceptance probabilities very 
close to unity) will in fact cost more than L2MC. Fortunately, this is irrelevant because the 
minimum cost for L2MC is far greater than the minimum for GHMC, the latter occurring 
for long trajectory lengths where our approximations are very good. 

Another implicit approximation we make is that we treat 6t and r as independent pa- 
rameters, although the trajectory length must be an integer multiple of the step size; again 
this is a very good approximation except for near to the Langevin limit. Of course we are 
also ignoring subtleties such as that a multistep leapfrog integration is cheaper than a series 
of single steps, that acceptance tests are a significant fraction of the cost for very short 
trajectories, and that HMC requires less memory than GHMC because the old momenta 
need not be saved. All of these facts would disfavour L2MC, but we shall see that it is not 
competitive even without these factors being taken into account. 



8.1 Linear Operators 

Applying equation (|3lD to equation (|40|) we find that the cost is minimised by choosing 
''^opt = 0, and at this value the Laplace transform of the single mode autocorrelation function 
becomes 

8.1.1 Exponentially Distributed Trajectory Lengths 

The integrated autocorrelation function for exponentially distributed trajectory lengths is 
obtained by substituting the results (|43|) for the integrals of equation (0) into Fm{P',^ = 

^I.e., the solution obtained by neglecting the dependence of P^cc on f when optimising the parameters d 
and f . 

''Setting — corresponds to never refreshing the momenta, and thus to a non-ergodic algorithm. This 
is just a peculiarity of linear operators in free field theory, and we can instead consider choosing § to be very 
small but non-zero to circumvent this difhculty. 
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■(^opt) and setting /3 = 0. The cost becomes 



m 6r 



(P_ - 2)^ ^ 4(P_ - 1) 



P., 



Pacc(Pacc-2)(^ 



Minimising this with respect to (f we find that 



opt 



1 - P. 



1 - iP 

2 ^cc 



and that the cost at the optimal parameters is 



41^(1 -Pa. 



opt J 



which corresponds to a dynamical critical exponent z = 1. 

For the lowest order leapfrog integration scheme where we must scale 5t oc V~'^ to keep 
Pace constant we thus find oc V^^ as expected. 



8.2 Quadratic Operators 

We can make the minimisation problem for quadratic operators manifestly algebraic by 
writing c = cos'i? in equation (^). The condition for the cost to be a minimum is that c^pt 
is a root of the polynomial 



^tl (^0,0 + ^0^2 + ^2^0) 



4 

opt 



2 Gi l '^0,2^2,0 + ^0,2 +^^0,2^^0,0) 



-'opt 



_L / O n+ 2^_ „ ^_(_ 3 r) ^+ 2^_ „ ^_|_ ^_|_ 2 „ ^_|_ 2^_|_ „ ^_(_ 2^_|_ ^ 2 



',2 y '^opt 

-G+2'c„p, + G+i' (49) 
lying in the interval [—1, 1]. 



8.2.1 Exponentially Distributed Trajectory Lengths 

Just as in equations ( pTD the extrema of the cost occur on the ideal defined by the polynomials 
d€./dc and d<L/dip. Finding the Grobner basis with respect to the purely lexicographical 
ordering with c < we find the point (Copt, v^opt) at which the cost is minimal is defined by 
the equations 

= (-4 + 3P_)c.p,^ + 4(l-2P_)(l-P_)c.p,3 

+ 16(1 - P_)c„p,2 + 4c„p, + P_ - 4 (50) 

= 4Pl(2-P.eJ(4-P_)(Pl-6P_ + 6)y..p,2 
+ (-4 + 3P..J(Pi + 2Pi - 10P..C + 8)c„p/ 
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+2(-l + P..J(4Pl + 6Pi - 47Pi + 68P_ - 32)c^J 
+ (-38Pl + 108Pi - 16P1 - 152P_ + 96)c„,,=^ 
+ (-16Pi + 56P1 - 20Pi - 168P1 + 248P_ - 96)c^,,' 
+ (35Pt - llOPi + 54Pl + 88P_ - 64)c.,, 

+8Pi - 60Pt + 126Pi - 62Pl - 48P_ + 32 (51) 



Using Sturm sequences we may easily show that equation ( pOf ) has exactly one real root in 
[0, 1] and none in [—1,0], and obviously equation (|5l|) has exactly one positive solution for 

V9„pt. The cost at the point (c„pt, v?opt) is given by 



V 



/ (7P_ - 3Pi - 4)^.,,'c^J + (-2P_ + l)c„p, + 1 

~l~ (2Pacc l)Copt Copt^ ~l~ ( Pace ~l~ ^)^opt 

V + (PL - 5Pacc + 4)Cop,¥5op,2 + (-4 + 3Pacc)Copt V. 



pt 

2 

opt 



P^^JtuiIP,,, - l)c„pt3</9„pt + (p,^,P,,jTm 



opt 

P ,o 

opt 



fopt 



This solution is a function of 5t and P^cc which are not independent variables, and using 
the results (0), and (^) we can compute the cost as a function of P^cc as shown in 
Figure ||. 

8.2.2 HMC 

The Hybrid Monte Carlo algorithm corresponds to setting d = tt/2, and thus we find that 



the optimal trajectory length in this case is (f^pt = l/v4 — Pa^c, corresponding to a cost 



2Vx/4 - P, 



This is also shown in Figure 0. 

8.3 Fixed Length Trajectories 

For fixed length trajectories we shall only analyse the case of L2MC for which the trajectory 
length (f = mdr. In this case Copt satisfies 



PaccCopt^(cOSV7„pJ^ + 3C„pt^Pacc + Pa,c (cOS V^ppt ) ^ - Pa,c - 2C„pt^ - 4c„pt Pace (cOS V9„pt ) ^ + 2 

Pace (c„pt2(cOSV?„pt)2 - c„p^2 _ (coS<^„pt)2 + 1) 



0. 



This, together with the corresponding cost, is also plotted in Figure ^. From this figure it 
is clear that the minimum cost occurs for P^^ very close to unity, where the scaling variable 
X = V5t^ is very small. We may then express Copt as 



X 



1/6 ^ ^(x\^/^ 11 3/a;\V2 



^-^'W +>IW -24^1f' + 
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Figure 2: Cost as a function of average Metropolis acceptance rate for the GHMC algo- 
rithm compared to HMC and L2MC for free field theory using the lowest order leapfrog 
integration scheme. The operator under consideration is the "magnetic susceptibility", i.e., 
the connected quadratic operator depending only on the lowest frequency mode. The cor- 
responding parameters, the momentum mixing angle i^^pt and the average trajectory length 
measured as a fraction of the correlation length ip^pt = Topt/'C ^-^e also shown, all as a function 
of the acceptance rate P^^^. The inset graph shows the region where the acceptance rate is 
very close to unity which is where the L2MC algorithm has its minimum cost. 
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and likewise 

P^^^ = 1 - J— 1 + + + 0(m^ 

V Stt V 20 800 ^ 

From these relations we find that the minimum cost for L2MC is 

/in \ 

Cp. = 2 (^) V'/'m-'/'. 

This result tells us that not only does the tuned L2MC algorithm have a dynamical critical 
exponent z = 3/2, but also it has a volume dependence of exactly the same form as HMC 
We may understand why this behaviour occurs rather than the naive V^/^m-i by the 
following simple argument. 

If Pace < 1 then the system will carry out a random walk backwards and forwards along a 
trajectory because the momentum, and thus the direction of travel, must be reversed upon a 
Metropolis rejection. A simple minded analysis is that the average time between rejections 
must be 0{l/m) in order to achieve z = 1. This time is approximately 



P.^rSr 1 



oo 

y 

n=0 acc '''' 



For small br we have 1 — Pj^c = erf \/kVbT^ oc yVbi^ where is a constant, and hence we 
must scale br so as to keep Vbr'^ jvr? fixed. Since the L2MC algorithm has a naive dynamical 
critical exponent -2=1, this means that the cost should vary as £ oc V(V"5r^m^^)^/^/m(5r = 



9 Autocorrelations and Frequency of Measurement 

In this final section we perform an elementary analysis of the general problem of determining 
the optimal frequency for making measurements of observables on a Markov chain. 

Suppose we wish to measure the expectation value (fi) of an operator by means of a 
Markov process (not necessarily HMC). If the costf] of making one Markov step is €5, and 
the cost of making one measurement of f2 is Ca/, how often should we make measurements 
in order to minimise the cost of the calculation? Due to the presence of correlations between 
successive measurements, the answer is not quite trivial. 

Consider a sequence of A^^ uncorrelated measurements of Q. The sample variance Vq is 
related to the intrinsic variance Vq in the distribution of Q in the usual way 

(n-M))^iL_pii,_a, (52) 



For the general case of Nc successive correlated measurements, from equation ([27|), the sample 
variance Vb is 

K = ((^ - {^)) ) = (1 + 2An) ^ ^ = (1 + 2An) ^ (53) 



*/.e., the cost measured in units of computer time. 
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so that on average 1 + 2Aq correlated measurements are needed to reduce the variance by 
the same amount as a single truly independent measurement. If the cost of measuring Q is 
small (large), then it should be beneficial to make more (less) than one measurement per 
decorrelation time. 

Let T be the total number of Markov steps required to generate one independent sample, 
and let K be the number of Markov steps between each measurement of Q, so that the 
number of measurements performed per independent sample is T/K. If we wish to make 
sufficient measurements to generate N independent samples, from equations (|5^) and ( ^31) 
we have 

where Ak is the integrated autocorrelation function for measurements separated by K 
Markov steps. The total cost €. of the entire process (steps plus measurements) is clearly 

(E = N(^€sT + (EM^y (55) 
Eliminating T from equations (^^ and (|55| ) we obtain 

€=N(^(ts + ^^ {1 + 2Ak)K. 

The integrated autocorrelation function Ak can be written in terms of the autocorrelation 
function Cq{KI) of section (|6.1|). Assuming equation ( |2^ ) we have 

oo CXD -| 

1=1 1=1 e I' ± 

and so 

€ = ms + ^m) [ ^k/n..,_, ) ■ 
Minimising with respect to K gives 

sinhx = X + a (56) 

where 

K , Cm 

X = — — and a ~ 



For small a where measurements are "cheap" we find 

1/3 



6€ 



M 



cxp . 



while for large a where they are "expensive", we obtain 

the slow logarithmic increase of -ft'opt being due to the exponential decay in the autocorrela- 
tion function Cn{l). The crossover point, K = 1, occurs when a ~ 0.1752. 
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10 Conclusions 



We have introduced some powerful techniques for analysing a wide class of algorithms for 
free field theory. Whereas inventing algorithms which are efficacious for free field theory 
but useless for interesting ones is a popular but fruitless exercise, the analysis of generic 
algorithms for free fields is apparently very informative. The reason that this is so is because 
our algorithms are sufficiently dumb that they spent most of their time in dealing with the 
trivial almost free "modes" of the system. 

As evidence of the more general applicability of our analysis we point to the excellent 
agreement with empirical Monte Carlo data of the erfc form for the Metropolis acceptance 



probability as a function of integration step size |2^. Futhermore preliminary results 
for simple models indicate that their autocorrelations have at least the same form as 
expected from our free field theory results. 

Perhaps the most surprising result is the failure of the L2MC (Kramers) algorithm to 
have superior performance to the HMC algorithm. Despite the fact that the noise is put 
into the Markov process more smoothly the desirable properties of the L2MC algorithm are 
upset by its "Zitterbewegung" caused by its rare Metropolis rejections. 

It is also somewhat unexpected that the HMC algorithm performs nearly as well as the 
full GHMC algorithm, of which it is a special case, when the parameters of both are chosen 
optimally. The broad minimum of the costs of these algorithms as a function of acceptance 
rate has the pleasant consequence that no very fine tuning of parameters is required. 

The result that the optimal HMC trajectory length is about the same as the correlation 
length of the underlying field theory is significant, in that it indicates that the temptation to 
use shorter, and hence cheaper, trajectories for systems near criticality should be avoided. 

The results pertaining to higher order (Campostrini) integrators p2| , p3| are of theoretical 
interest, but in practice they do not seem to be very useful because the trajectories for 
interacting field theories are chaotic p7| , p8| and apparently limited by the intrinsic instability 
of the integrators [O]. 
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A Inverse Laplace Transforms 
for Autocorrelation Functions 



A.l Exponentially-Distributed Trajectory Lengths 

In this case the Laplace transform of the autocorrelation function, F{P) — CC{t), is a rational 
function with the numerator of lower degree than the denominator. If the denominator is 
square free then we have the partial fraction expansion 




Since 

£6""* = / dt e-^*e-"* = (Re a + /3 > 0) , 

Jo a-\- p 

we have 

k 

The roots of the denominator of F come in complex conjugate pairs because the coeffi- 
cients in F are real. The autocorrelation function can therefore always be written as a sum 
of real exponentials (corresponding to the real roots) and of real exponentials multiplied by 
cosines (corresponding to pairs of complex conjugate roots). 

If the denominator of F is not square free then the repeated roots give terms of the form 
A/ [13 — BY in the partial fraction expansion. Consider then 

/■oo 

£^n-lg-at = / ^^g-/3t^n-lg-at ^ (q; + /3)~"r(n) (Rc a + /3 > 0), 

Jo 

This serves to give the inverse Laplace transform in the general case. 



A. 2 Fixed Length Trajectories 

Here F is a rational function of C = e"'^^ with the numerator being of lower or equal degree 
than the denominator. For simplicity we first consider the case where the denominator is 
square free, whence by partial fractions 

F = C + ^ ^' 



k C-Bk 

Observe that _ 

/■oo 

£5(t) = / e-^'5(t) = 1, 
Jo 

and, more generally, 

oo oo -1 

j=0 j=0 ^ 
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so we have 

A, °° 

C'^F = C6{t)-j:-^j:B,^6{t-rj) 



j=0 



k 3=0 

All the the roots must satisfy > 1 for the geometric series to converge at /3 = (C = 1), 
which is necessary for the integrated autocorrelation function -F(O) to be finite. 

In the general case where the denominator of F has multiple roots the general form of 
the inverse Laplace transform may be obtained from the identity 

j=Q \ J / 

A. 3 Example of Computation of Autocorrelation Function 

To illustrate the computation of autocorrelation functions we shall explicitly evaluate the 
connected autocorrelation function for for the HMC algorithm with exponentially dis- 
tributed trajectory lengths. The Laplace transform of the autocorrelation function is 

/33 + 2r/?2 + (r2 + 4m2)/3 + 2mV 
If we write this in terms of two distinct roots /3i and P2 of the denominator 

r (r2 + 2r/3 + /32 + 2m2) 

we can expand this in terms of partial fractions to give 

r (2 + 4 r^Pi + (2 (3i^ - m^) + p^r + 32 + 4 /^i^m^) 
^^^^ " (/3 - (2 - 13 m^r^ + 64 m^) 

-4 r^(3i + (-2 (3i^ - 13 m2 - 2 (32P1) 
+ (-7 /32 - 8 m^r + 16 + (-4 - 4 
((3 - (32) (2 - 13 m^r^ + 64 m^) 
r ((2 /32/3i + m^) + (7/32 + 7 A) m^r + 16 + 4 (32(3im'^) 
(2r + (3i + (32 + (3) (2H- 13mV2 + 64m4) ' 
The inverse Laplace transform of this is immediately obvious 

r (2 + 4 r^(3i + (2 - m^) + /3ir + 32 + 4 /^i^m^) e^i* 



2 

m 



+ 



+ 



+ 



+ 



2 - 13 mV2 + 64 
-4 r^A + (-2 - 13 m2 - 2 /32/3i) ^^^^ 
' +(-7/32-8/3i)mV + 16m^+(-4/3i^-4/32/3i)m2 ^ 
2H - 13 m^r^ + 64 

(2/32/3i + m2) r2 + {7(32 + 7 A) m^r >| ^_(2.+/3i+/3.)t 

+ 16 m'^ + 4/32/3im/ 



2r4 - 13 m2r2 + 64 
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Figure 3: The autocorrelation function and the cumulative cost are shown as functions of 
MD time for the optimal trajectory length and those costing 50% more. 

In figure ^ we show this function for m = 10~^ and the optimal value r^p^ = \/3m, together 
with the two neighbouring values of r = {3^/3± vT5)m/2 for which the cost is 50% greater. 

B Evaluation of Spectral Averages 

In this section we show how to evaluate the following spectral sums for one-dimensional free 
field theory: 
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a, 



= (59) 



,2 



Cm 2 



IB')"l'-."(.!y|. .« 

i ^ cos(2Ta;p) (61) 

E^4^ H^t/ cir,^i EW4r), (62) 
13=0 P- V7=o / fceZ 

i ^ <sin2(ra;,) (63) 

1 1 / 1 r . 



where 



^J = m2 + 4sinM^). (65) 



B.l Completeness relations 



Let ^ be a group with a family of D-dimensional representations (J)g : Q ^ (C ^ C ), z.e., 
0g(5') : C'^ C^, labeUed by some parameter q; thus (f)q{gh) = 4>q{g)(f)q{h) for g,h ^ Q. 
Summing over group elements (or integrating with respect to Haar measure in the continuous 
case) gives 



(Tq 



g 



and multiplying by (f)q{h) for any h & Q we find 

(f>q{h)aq = (f)q{h) / d/i (f)q{g) = / dn(t)q{h)(pq{g) 

JQ JQ 

^ dfx (f)q{hg) ^ dfx (j)q{g') = aq 

JQ JQ 

where we have made use of the (left) invariance of Haar measure d^. Hence either 0g(/i) = 
1 V/i G ^ or (jg = 0. If we let g = label the identity representation of and all other (pq 
be non-identity representations, we have 

Gq = 8{q). 

Example B.1.1 Let Q ^Zy, qe ly, and (t)q{p) = e^'^^P^/^; 

p&Lv 
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Example B.1.2 Letg^Z, qeSi^ [0, 27r], and (j)g{p) = e'Pi: 
Example B.1.3 Let G ^E>i, qeZ, and (j)g{p) = e'Pi: 



B.2 A Poisson Resummation Formula 

Let / : §1 ^ Si be a periodic function. As it is periodic it obviously wants to be expanded 
in eigenfunctions of §i: 

27r 

IJ 



fip) = / dp'S{p'-p)f{p') [pes 
Jo 

p2'K 1 

= I dp'^Y.^'^'^'-'^m 



1 f'^'^ 



t^Z 27r Jo 



hence the "spectral average" 



As / is real valued we may take the real part of this identity to obtain 

<f) = 7 E / i^-y) = E ^ fj dp' cos{p'kV)m. (66) 
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B.3 Multiple Angle Expansion 

The following simple identity is often useful 



sin^" e 



2i V 

1 \ 2" ^ f2a 
2i 



2a 



{a E N) 



'_)9e2j(a-g) 



g=o \ y / 

^ ' ^ ^^"^ i-r + E (^f ) (-)^2 cos[2(a - g)^] 



2i/ \ \ a 
1 f /2a 



2a 



22c 



a 



+ E L .J(-)'''2cos(2g'^ 



g'=l 



, a — g 



(67) 



B.4 Free-field Spectral Sums 

Let us first consider equation (|57D . From equation ( |66D we have 



^ E cos [^Y^j = E ^/ cos(p'fc\/) cos(pV) 
1 /-s^ 1 



E 



27r 7o 2 



|cos[]9'(A;V^ + q')] + cos[p{kV — 



fceZ 



}■ 



Using the multiple-angle expansion ([67| ) we obtain 



0" 



(a) 



E 



sm 



peZv 

- E - 

V 22° 



2a 

V 
'2a 



a 



1 

22° 

1 

22° 

1 

22^ 



^2a^ 
^2a^ 
^2a^ 



2a 



.'1 



E cos —— 



peZv 



+ E 



2a 



E (-)'^ 

k=l 



' 2a ^ 
,a — kV I 



' 2a ^ 
^a + A;\/^ 



From this result we may obtain an expression for the quantity a^J, defined in equation 



m + 4 sm I — 
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/3=0 



1(7) 



Example B.4.1 



an 



^ E = (6 + + m') + 6v,i (-6 - 4m') + 6v,2 ■ 2; 



Example B.4.2 



(m^ + 6m^ + ISm^ + 20) 
(-6m^ - ISm^ - 20) + 5^,2 (em^ + 12) - 5y,3 ■ 2. 



- y ^6 



In order to evaluate Equation (|6T|) we first investigate the simpler case where m = 
^0 = ^ E cos(2™p) 



1 /■^'^ 

— / dp cos(p'kV) cos 
keZ 2^ ^0 



27r 7o 
1 r 



At sin 



y - r dO cos(2A;r^)cos[4rsin^] 
71 Jo 



keZ 
E J2kv(4r). 



For small m we Taylor expand 



m2=0 



Consider the quantity 



1 d d a'm2 

2 (im^ c/r r 



1 / 1 1 



2 c/r {rdrn'^V 



cos(2ra;p) 



d ( I ^ , ,1 

- - E --(2-p)^ 

^ y cos{2TUJp) = a^2; 



46 



upon integrating with respect to r we obtain 



dm? 



-2t I dr' am?, 



and hence 



d" 



d{rr?)P 



m?=0 



(-2)^r / dri n dT2 . . . r^-i / dr/^ao 



By this means we have obtained the desired result 



"^mZ = 2^ -T, Ml '^7 / dT^+i (Jo. 



(3=0 f^^- \7=0 

Equation (|63D may be found by differentiating with respect to r, 



V 



Example B.4.3 



^ E ^f{l-cos(2™,)} 
1 H 1 / 1 rfMV 



an 



— J2 ^pSin2(ru;p) 



E 

/3=0 



/5! V 4rfrV d(m2)/5 



m^=0 



2 2/3! 



fc6Z/3=0 



X 



^(3-1 



X I n ^7 / ^^7+1 I J2kv(4r). 

17=0 



The f3 = term in the sum is 



feeZ \ ' 
From the recurrence relation 



-\ E f-l:^V J2w(4r) = -8 E W'"(4r). 



fceZ 



Jn — "(Jn-l ~ Jn+1, 



we find that 



''Jn-k+ 27(^)5 
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which we may prove by induction: 



1 ^ (^\( ^7 '^Jn-k+27(^ 
Ok 2-^ 



7=0 



dz 



1 ^ (k\ 

2fcTT (-)^(jn-k+27-l(^) - Jn-k+27+l(^)) 

Yl L (-)'^(jn-(k+l)+27(^) - Jn-(k+l)+2(7+l)(2:)) 
7=0 \ // 

-)^Jn-(k+l)+27(2;) + 



1 A (k^ 



2k+l 



2k+ 



k 

HE 

7=1 



' k ^ 
.7-1. 



1 + 1 



2fc+i 



+Jn-(k=l)(2;) — Jn+(k+l)(2:) 
(-)'^Jn-(k+l)+27(2:), 



7=0 



7 



using (^) + (^*: J = (''+^) . l^e therefore find that 

-8 E w"(4r) = -2 ^ H (-rj 



2kV-4+27 



(4r) 



fceZ fceZ7= 

-2 E {J2kv-4(4r) - 4J2kv-2(4r) + 

feeZ 

+6J2kv(4r) - 4J2kv+2(4r) + J2kv+4(4r)} 
-3Jo(4r)+4J2(4r)- J4(4r)- 

oo 

-2 E{j2kv-4(4r) - 4J2kv-2(4r) + 



k=l 



+6J2kv(4r) - 4J2kv+2(4r) + J2kv+4(4r)|. 



The (3=1 term in the sum is 



H ^ ■ ^ y J2kv(^') (where f = 4r) 

fceZ 



XI 4^ I / J2kv(T"') + 

m2E{4J2kv"(r) + fJ2kv'"(r)} 

fceZ 



2kVlT") 



EE 



)'^J2kV-2+27('r) + 
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3 /gN 



= X]{j2kV-2(T) - 2J2kv(T) + J2kV+2(T) + 

3f, ... 3f 



feeZ 

+ gJ2kV-3(T) - yJ2kV-l(T) + yJ2kV+l(T) - ^J2kV+3(T)| 



= m2|-2Jo(4T) + 2J2(4t) + 3tJi(4t) - tJ3(4t) + 

oo 

+ E(2J2kv-2(4r) - 4J2kv(4r) + 2J2kv+2(4r) + 

k=l 

+rJ2kv-3(4r) - 3TJ2kv-i(4T) + 
+3TJ2kv+i(4T) - rJ2kv+3(4r))|. 

B.5 Results for Campostrini 

35 - 35Jo(4t) + 56J2(4t) - 28J4(4t) + 8J6(4t) - J8(4r 
40 - 40Jo(4t) + 64tJi(4T) + 31J2(4r) - 8J4(4r) + J6(4t) 













+ 




+ 


an 










+ 




+ 



m 



18 + (IGi^ - 18)Jo(4t) + 48iJi(4T) + 8J2(4t) - J4(4T)]m^ + O(m^), 

462 - 462 Jo (4r ) + 792 J2 (4r ) - 495 J4 (4r ) 
+220J6(4r) - 66J8(4r) + 12Jio(4r) - Ji2(4r 

756 - 756Jo(4r) + 1024tJi(4r) + 698J2(4r) 
-256J4(4r) + 69J6(4r) - 12JsW) + Jio(4r) 

525 + (256^2 - 525) Jo(4t) + 1152Ui(4r) 
+316J2(4t) - 74J4(4t) + 12J6(4t) - J8(4t) 



m 



+ 



+ 



6435 - 6435Jo(4r) + 11440J2(4r) - 8008J4(4r 
+4368J6(4t) - 1820J8(4r) + 560Jio(4r) 
-120Ji2(4t) + 16Ji4(4t) - Ji6(4t) 

13728 - 13728Jo(4t) + 16384Ui(4t) 

+ 14075J2(4r) - 6128J4(4r) + 2185J6(4r) 
-608J8(4r) + 123Jio(4r) - 16Ji2(4r) + Ji4(4t) _ 

12936 + (4096t2 - 12936)Jo(4r) 
+24576tJi(4r) + 9280J2(4r) - 2807J4(4r) 
+688J6(4r) - 128J8(4r) + 16Jio(4r) - Ji2(4r 

92378 - 92378Jo(4r) + 167960J2(4r) 
-125970J4(4r) + 77520J6(4r) - 38760J8(4r) 
+15504Jio(4t) - 4845Ji2(4t) + 1140Ji4(4t) 
-190Ji6(4t) + 20Ji8(4t) - J2o(4r) 



m 
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+ 



+ 



243100 - 243100Jo(4r) + 262144tJi(4r) 
+267814J2(4r) - 129872J4(4r) + 54252J6(4r) 
-19024J8(4r) + 5420Jio(4r) - 1200Ji2(4r) 
+ 193Ji4(4r) - 20Ji6(4r) + Ji8(4r) 

289575 + (65536^2 - 289575)Jo(4r) 
+491520tJi(4r) + 233612 J2(4t) - 82714J4(4r) 
+25164J6(4r) - 6392J8(4r) + 1300Jio(4r) 
-198Ji2(4t) + 20Ji4(4t) - Ji6(4t) 



m 



C Spectral sum for one dimensional free field theory 

We wish to evaluate the spectral sum 



V 



peZv 



for large V . Using the Poisson resummation formula we obtain 
Upon substituting z = e*^' we get 



rm? + 4sin^(ip') 



^kV + ^-kV 



1^ 47r i|=i z^ - {m? + 2)z + 1 



E 



dz 



27r /|z|=i {z - Hj^){z - ^x-) m\/nn? + 4 



V 



l + /i 



where = 1 + ± \m\/m? + 4. For m <C 1 wc have /i_ = 1 — m + 0{w?), and thus 



coth imV 



rn^Jm? + 4 

which is correct to all orders in the large V asymptotic expansion. 



D Spectral sum for two dimensional free field theory 



We want to evaluate the spectral sum 



E 



where 



uj^ = m + 4 sm — h sm 



2Pj/7r 
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Using Poisson resummation we find that for L^, Ly ^ oo 



(27r)2 Jo 



2-K 



+ 4 ( sin^ — + sin^ — 



so upon substituting z = e*^^ we obtain 



a — 



47r2 
i 



2^ dz 

kl=i z^ - (m2 + 4sin2(ipp)z + 1 

<- dz 



dp'y 

dp'y 



An^ Jo ^ J\z\=i {z - ijl){z - l/fji) 



where 



u = — + 2 sin^ ^ + 1 
^2 2 



Evaluating the contour integral gives 



1 

dp'y 



jj, J Att Jo 

We now make the substitution z' = e^^'y , which leads to 

i f dz' 



m 

T 



+ 2sin2?^ + ll -1 



cr 



27r y|z'|=i ^^(^^2 + 4)^/ _ _ if _ 4^/2 
i r dz' 



27r J\z'\=i J^z' - a){z' - l/a){z' - b){z' - 1/6) 



where 



= ^ + 3-^(^ + 3) -l = (3-2x/2)-^(3x/2-4)m2 + 



^ + Ij - 1 = 1 - m + + 



Since l/a>l/6>l>6>a>0 the contour integral may be shrunk to be the integral 
along the branch cut from a to b, 



1 

a — — 



dz' 



^J{z' -a){b- z'){l/a - z'){l/b - z') ' 



If we change variables to 



X = 



(a + l)(;g^-l) 
(a-l)(^' + l) 
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we obtain 

a = 



7T{a + l){b-l) Ji y^(a;2 - 1)(1 - k'^x^) 



where k' = . With the substitution ^ = jJl — ^ we find that 



(a+l)(6-l)- 



^(a;2 - 1)(1 - A;'2a;2) -/o -^2)(i _ p^2) 

where is the complete elliptic integral and /c^ + /c'^ = 1; thus 

2v^ 



K(A;) 



a 



-K{k). 



7r(a + 1)' 
This may be further simphfied using Landen's transformation 



with r] = which leads to 



( b — a 

-K 



7r(l — ah) \1 — ah ) 
using the fact that K is an even function. Expanding in powers of m we obtain 

a = -(l + 0{m))K(l - V2m + Oim^)) = — Inf^) +0(mlnm), 



because 

k 



Kik) = K(Vl - k'^) ~ In ^ as fc' ^ 0. 
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